"""
GM1050 Advanced Principles and Algorithm of GIS
2010/2011

- Mini Project - The simplification of a map

A shapefile of departments in western Europe is given.
By the Douglas-Peuker algorithm this shapefile is simplified.
Where several topological relations should remain

@author: Bas, Josafat and Elise
"""
#-- general import
import math
import os
#from Vertex import Vertex
from dcel1 import DCEL
#from Douglas_Peucker import algorithm
from Douglas_Peucker2 import simplify_points
#-- OGR 
try:
    from osgeo import ogr
except ImportError:
    import ogr
#-- Shapely
from shapely.wkb import loads as wkbloads
from shapely.wkt import loads as wktloads

from shapely.geometry import LineString





     

def main():

    createPolygons();

def createPolygons():
      #1) Load all polygons
    #shapefile = 'final.shp';
    #shapefile = 'departementsmall.shp';
    shapefile = 'europe_nuts.shp';
    #shapefile = 'europe_nuts_mini2.shp';

    print "reading: " + shapefile
    
    shapelyFeatures = read_shapefile(shapefile)

    print "Number of Features: " + str(len(shapelyFeatures))

    polygons = simplifyPolygons(shapelyFeatures)# i just recieve a list of indexed coordinates
    
    wktpolygons = createPolygonsWkt(polygons)

    outfile = raw_input('give filename (+.shp) of output file ... ')
    writePolygonShapefile(shapelyFeatures, wktpolygons, outfile)
    #construct_graph(polygons)
  
def createLines():
        #1) Load all polygons
    #shapefile = 'final.shp';
    #shapefile = 'departementsmall.shp';
    shapefile = 'europe_nuts.shp';
    #shapefile = 'europe_nuts_mini2.shp';
    print "reading: " + shapefile
    shapelyPolygons = read_shapefile(shapefile)
    print "Number of Features: " + str(len(shapelyPolygons))
    
    #2) Extract edges and nodes - (not line segments, but a graph).
    
    #3) Structure for each polygon which edges is bounding him

    polygons = simplifyPolygons(shapelyPolygons)
    
    text = createPolygonsWkt(polygons)

    outfile = raw_input('give filename (+.shp) of output file ... ')
    write_shapefile(text, outfile)
    #construct_graph(polygons)
    

def read_shapefile(infile):
    """
    Function that reads a shapefile and prints some info about it.
    
    Input:
        infile:   shapefile to read
    Output:
        - A dicctionary

    """
    print "Reading", infile 
    ge = ogr.Open(infile)
    layer = ge.GetLayer(0) #-- a simple shapefile always has 1 layer only
    ld = layer.GetLayerDefn()
    
    #-- Find and print the geometry type of the shapefile
    print ld.GetGeomType()
    s = "Geometry type: "
    if ld.GetGeomType() == ogr.wkbPoint:
        s += "Point"
    elif ld.GetGeomType() == ogr.wkbLineString:
        s += "LineString"
    elif ld.GetGeomType() == ogr.wkbPolygon:
        s += "Polygon"
    else:
        s += "Unknown"
    print s
    
    #-- Find and print the number of geometries
    noFeat = layer.GetFeatureCount()
    print "Number of features:", noFeat
        
    #-- Find and print the attributes
    print "Attributes:"
    
    
    
    listOfAttributeTypes = []
    
    for i in range(ld.GetFieldCount()):
        listOfAttributeTypes.append(ld.GetFieldDefn(i).GetNameRef())#add the attribute name to the list
        print "\t", ld.GetFieldDefn(i).GetNameRef() + " " 
    print "The features are:"
   
    lsFeatures = [] #-- create an empty list of Features (of geometries)
    for i in xrange(noFeat):

        f = layer.GetFeature(i)
        wkt = f.GetGeometryRef().ExportToWkt() #-- Well-Known Text

        #print str(i) + "\t" + wkt 

        
        pshape = wktloads(wkt)

        feature = {'id': str(i), 'geom': pshape  , 'fieldNames' :listOfAttributeTypes}  
        for j in xrange(len(listOfAttributeTypes)):
            feature[ str(j) ] = f.GetFieldAsString(j); 
        
        lsFeatures.append(feature)
    
    return lsFeatures

#2) Extract edges and nodes - (not line segments, but a graph).

def simplifyPolygons(shapelyPolygons):

    dcel = DCEL()
    dcel.convert(shapelyPolygons)## just import the simple features into DCEL
    
    lineSegments = dcel.findLineSegments();## the line segments are a list of tuples
    
    tolerance = float(raw_input('Please give tolerance (e.g.: 0.03) [deg] '))
    for lineSegmentKey in iter(lineSegments):
        lineSegment = lineSegments[lineSegmentKey]
        oldline = dcel.convertLine2Tuple(lineSegment);
        newline = simplify_points(oldline, tolerance);
        dcel.replaceLine(oldline, newline);############check for errors

    # this call will return a dictionary
    # {face = ''}
    polygons = dcel.createPolygons2()
    
    ## converting  a dict 'polyid' , [vertexes] to 'polyid'


    return polygons


def simplifyLines(shapelyPolygons):

    dcel = DCEL()
    
    dcel.convert(shapelyPolygons)

    lineSegments = dcel.findLineSegments();
    
    
    tolerance = float(raw_input('Please give percentage to keep ... [0...1]'))
    #percentage = float(raw_input('Give the percentage of points to be kept [0...1] ...'))
    lines = []
    for lineSegmentKey in iter(lineSegments):
        lineSegment = lineSegments[lineSegmentKey]
        #print "Segment: " + lineSegmentKey 
        oldline = dcel.convertLine2Tuple(lineSegment);

        newline = simplify_points(oldline, tolerance);

        dcel.replaceLine(oldline, newline);############check for errors
        lines.append(newline)
    

    return lines

def createLinesWkt(lines):
    wkt = []
    for i in range(len(lines)):
        # for each polyline
        line = ''
        for j in range(len(lines[i])):
            # for each coordinate
            line = line + str(lines[i][j][0])\
            + ' ' + str(lines[i][j][1]) + ' , '
        polyline = 'LINESTRING (' + line[0:-2] + ')'
        
        wkt.append(polyline)
    return wkt




def createPolygonsWkt(polygons):
    ''' 
        I should recieve a dictionary with polygons
    '''
    wkt = []

    #print polygons
    for i in xrange(len(polygons)):

        polyline = 'POLYGON ((' 
        polygon = polygons[str(i)]
        tuple1 = polygon.pop()
        #print tuple1
        polyline = polyline + str(tuple1[0]) + " " + str(tuple1[1])
        
        for tuple in polygon:
            polyline = polyline + " , " + str(tuple[0]) + " " + str(tuple[1])
        
        polyline = polyline + "))"
        wkt.append(polyline)
    return wkt

def writeLineShapefile(lsFeat, outfile):
    """
    Create a new shapefile and writes it to the disk.
    
    Input:
        lsFeat:    list of Features (a dico with geom/name/pop)
        outfile:   path for the new file to create
    Output:
        - (the shapefile is written to the disk)

    """
    driver = ogr.GetDriverByName('ESRI Shapefile') #-- we create a new SHP file
    if os.path.exists(outfile):
            driver.DeleteDataSource(outfile) #-- if it exists, overwrite it
    ds = driver.CreateDataSource(outfile)
    layer = ds.CreateLayer(outfile, geom_type = ogr.wkbLineString) #-- we create a SHP with polygons
    
    #-- create the 1 attributes (called Fields in OGR jargon) "name" and "population"
    fd = ogr.FieldDefn()
    layer.CreateField(fd)
        
    for i in lsFeat:
        f = ogr.Feature(feature_def = layer.GetLayerDefn()) 
        p = ogr.CreateGeometryFromWkt(i)
        f.SetGeometry(p)
        layer.CreateFeature(f)
        f.Destroy()
    ds.Destroy()
    print "\nShapefile saved as:", outfile

def writePolygonShapefile(shapelyFeatures, wktpolygons , outfile):
    """
    Create a new shapefile and writes it to the disk.
    
    Input:
        lsFeat:    list of Features (a dico with geom/name/pop)
        outfile:   path for the new file to create
    Output:
        - (the shapefile is written to the disk)

    """
    driver = ogr.GetDriverByName('ESRI Shapefile') #-- we create a new SHP file
    if os.path.exists(outfile):
            driver.DeleteDataSource(outfile) #-- if it exists, overwrite it
    ds = driver.CreateDataSource(outfile)
    layer = ds.CreateLayer(outfile, geom_type = ogr.wkbPolygon) #-- we create a SHP with polygons
    
    #-- create the 1 attributes (called Fields in OGR jargon) "name" and "population"

    
    fd = ogr.FieldDefn()
    
    feat = shapelyFeatures[0]
    names = feat['fieldNames']
    print "-----------"
    print names
    print "-----------"
    for i in names:
        print i
        fd.SetName(i)
        fd.SetType(ogr.OFTString)  #-- type is String
        layer.CreateField(fd)
        
    
        
    for i in xrange(len(wktpolygons)):
        f = ogr.Feature(feature_def = layer.GetLayerDefn()) 
        p = ogr.CreateGeometryFromWkt(wktpolygons[i])
        
        f.SetGeometry(p)
        feature = shapelyFeatures[i]
        attributeNames = feature['fieldNames']

        for j in xrange(len(attributeNames)):
            f.SetField(j , feature[str(j)]) 
        

        layer.CreateFeature(f)
        f.Destroy()
    ds.Destroy()
    print "\nShapefile saved as:", outfile







if __name__ == "__main__":
    print "Running"
    main()
    #raw_input('press enter')


